#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>
#include<time.h>
#include<direct.h>
#include<sys\\stat.h>
#include<sys\\types.h>

#define mSol 1
#define mTerra 3.003167261575587e-6
#define mJupiter 9.543512141169373e-4

#define PI 3.141592653589793
//#define G 2.959572553206070e-004
#define G 4*PI*PI
#define diaAno 365.25
//#define G 1
#define perielioTerra 0.983235919814466
#define vPerielioTerra 6.389661159796174
#define excentricidadeTerra 0.0167

void exteriores(void);
void inicializacao(void);
void condicoesIniciais(void);
void coordenadasCentroMassa(void);
void criarArquivos(void);
void escreverArquivos(void);
void ruth(void);
void calcularDistancia(double *v,double *valorNorma, int dimensao, int instante, int T,int corpos);
void Fr(double *dx, double *p, double *m, int inicio,int dimensao, int T);
void ruthR(int dimensao, double *r,double dt,double *dx,double c,int inicio,int T);
void Fp(double *dx,double *r, double *m,double *distancia, int inicio, int dimensao,int T, int corpos);
void ruthP(int dimensao, double *p,double dt,double *dx,double d, int inicio, int T);
void calcularTempo(void);
void fecharArquivos(void);

double  e,a2,*m, *r, *p,   tFinal,    dt,     c[4]  ,  d[4] , expoente, temp,*a,*b,*c2,*dx, *distancia,tempo=0;
int dimensao, n, T, corpos,intervaloGravacao,minutos=0,segundos=0,horas=0,tempo2=0;
//2*N+1
FILE *f[13];
char *exemplo,*aux,*SO,*tipoArquivo;
clock_t inicio,fim,total;

int main(int argc, char *argv[])
{
	inicio=clock();

	if(argc==7)
	{
		SO=argv[1];
		intervaloGravacao=atoi(argv[2]);
		e=atof(argv[3]);
		a2=atof(argv[4]);
		tFinal=atof(argv[5]);
		n=atof(argv[6]);
	}
	else
	{
		SO=argv[1];
		intervaloGravacao=atoi(argv[2]);
		tFinal=atof(argv[3]);
		n=atof(argv[4]);
	}

	inicializacao();

	criarArquivos();

	escreverArquivos();

	ruth();

	calcularTempo();

	fecharArquivos();

	return 0;
}

void inicializacao(void)
{
	dimensao=3;
	corpos=6;
	T=(int)pow(2.0,n)+1;
	dt=tFinal/(T-1);
	printf("%lf\n\n",dt);

	condicoesIniciais();

	//coeficientes do metodo
	expoente=1.0/3;
	temp=2-pow(2,expoente);
	temp=1/temp;

	c[0]=temp/2;
	c[3]=c[0];
	c[1]=temp/2;
	c[1]=c[1]*(1-pow(2,expoente));
	c[2]=c[1];
	d[0]=temp;
	d[2]=temp;
	d[1]=temp;
	d[1]=-1*d[1]*pow(2,expoente);
	d[3]=0;
}

void condicoesIniciais(void)
{
	int i;

	m =(double *)calloc(corpos, sizeof(double));
	r =(double *)calloc(corpos*dimensao, sizeof(double));
	p =(double *)calloc(corpos*dimensao, sizeof(double));
	dx =(double *)calloc(dimensao*corpos, sizeof(double));
	distancia =(double *)calloc(corpos*corpos, sizeof(double));

	exteriores();
	coordenadasCentroMassa();

	//velocidade centro de massa ser nula
	for(i=0;i<dimensao;i++)
	{
		p[0*dimensao+i]=-p[1*dimensao+i]-p[2*dimensao+i]-p[3*dimensao+i]-p[4*dimensao+i]-p[5*dimensao+i];
	}
}

void coordenadasCentroMassa(void)
{
	int i,j;

	double M=0;
	double *rcm;
	rcm =(double *)calloc(dimensao, sizeof(double));

	for(i=0;i<corpos;i++)
	{
		M=M+m[i];
	}

	for(i=0;i<dimensao;i++)
	{
		rcm[i]=0;
		for(j=0;j<corpos;j++)
		{
			rcm[i]=rcm[i]+r[j*dimensao+i]*m[j];
		}
		rcm[i]=rcm[i]/M;
	}

	//transformacao
	for(i=0;i<corpos;i++)
	{
		for(j=0;j<dimensao;j++)
		{
			r[i*dimensao+j]= r[i*dimensao+j]-rcm[j];
		}
	}
}

void exteriores(void)
{
	exemplo="C:\\dados";
	mkdir(exemplo);
	exemplo="C:\\dados\\JSUNP";
	mkdir(exemplo);

	//Sol
	m[0]=1.00000597682;

	//Jupiter
	m[1]=0.000954786104043;
	r[1*dimensao+0]=-3.5023653;
	r[1*dimensao+1]=-3.8169847;
	r[1*dimensao+2]=-1.5507963;
	p[1*dimensao+0]=0.00565429*m[1]*diaAno;
	p[1*dimensao+1]=-0.00412490*m[1]*diaAno;
	p[1*dimensao+2]=-0.00190589*m[1]*diaAno;

	//Saturno
	m[2]=0.000285583733151;
	r[2*dimensao+0]=9.0755314;
	r[2*dimensao+1]=-3.0458353;
	r[2*dimensao+2]=-1.6483708;
	p[2*dimensao+0]=0.00168318*m[2]*diaAno;
	p[2*dimensao+1]=0.00483525*m[2]*diaAno;
	p[2*dimensao+2]=0.00192462*m[2]*diaAno;

	//Urano
	m[3]=0.0000437273164546;
	r[3*dimensao+0]=8.3101420;
	r[3*dimensao+1]=-16.2901086;
	r[3*dimensao+2]=-7.2521278;
	p[3*dimensao+0]=0.00354178*m[3]*diaAno;
	p[3*dimensao+1]=0.00137102*m[3]*diaAno;
	p[3*dimensao+2]=0.00055029*m[3]*diaAno;

	//Netuno
	m[4]=0.0000517759138449;
	r[4*dimensao+0]=11.4707666;
	r[4*dimensao+1]=-25.7294829;
	r[4*dimensao+2]=-10.8169456;
	p[4*dimensao+0]=0.00288930*m[4]*diaAno;
	p[4*dimensao+1]=0.00114527*m[4]*diaAno;
	p[4*dimensao+2]=0.00039677*m[4]*diaAno;

	//Plutao
	m[5]=1.0/(1.3e8);
	r[5*dimensao+0]=-15.5387357;
	r[5*dimensao+1]=-25.2225594;
	r[5*dimensao+2]=-3.1902382;
	p[5*dimensao+0]=0.00276725*m[5]*diaAno;
	p[5*dimensao+1]=-0.00170702*m[5]*diaAno;
	p[5*dimensao+2]=-0.00136504*m[5]*diaAno;

}

void criarArquivos(void)
{
	char *raiz,*diretorio,*aux,*barra;
	int i;

	raiz=(char*)calloc(100,sizeof(char));
	aux=(char*)calloc(100,sizeof(char));
	diretorio=(char*)calloc(100,sizeof(char));

	if(!strcmp(SO,"w"))
	{
		barra="\\";
	}
	if(!strcmp(SO,"l"))
	{
		barra="/";
	}

	strcpy(raiz,exemplo);
	strcat(raiz,barra);
	//mkdir(raiz,S_IRWXU);
	mkdir(raiz);
	sprintf(aux,"%0.0lf",tFinal);
	strcat(raiz,aux);
	sprintf(aux,"%d",n);
	strcat(raiz,barra);
	//mkdir(raiz,S_IRWXU);
	mkdir(raiz);
	strcat(raiz,aux);
	strcat(raiz,barra);
	//mkdir(raiz,S_IRWXU);
	mkdir(raiz);

	for(i=0;i<corpos;i++)
	{
		strcpy(diretorio,raiz);
		strcat(diretorio,"r");
		sprintf(aux,"%d",(i+1));
		strcat(diretorio,aux);
		strcat(diretorio,".txt");
		f[i] = fopen(diretorio,"w");
		if(f[i] == NULL)
		{
			printf("Erro ao abrir arquivo!!!\n");
			exit(1);
		}
		strcpy(diretorio,raiz);
		strcat(diretorio,"p");
		sprintf(aux,"%d",(i+1));
		strcat(diretorio,aux);
		strcat(diretorio,".txt");
		f[i+corpos] = fopen(diretorio,"w");
		if(f[i+corpos] == NULL)
		{
			printf("Erro ao abrir arquivo!!!\n");
			exit(1);
		}
	}

	strcpy(diretorio,raiz);
	strcat(diretorio,"tempo");
	strcat(diretorio,".txt");
	f[2*corpos] = fopen(diretorio,"w");
	if(f[2*corpos] == NULL)
	{
		printf("Erro ao abrir arquivo5!!!\n");
		exit(1);
	}


	free(aux);
	free(raiz);
	free(diretorio);
}

void ruth(void)
{
	int i,j;

	for(i=1;i<T;i++)
	{
		//printf("%lf%\n",(double)i/(T-1));
		//copiar(r,p,i,dimensao,T,corpos);
		for(j=0;j<4;j++)
		{
			Fr(dx,p,m,i*dimensao,dimensao,T);
			ruthR(dimensao,r,dt,dx,c[j],i*dimensao,T);
			calcularDistancia(r,distancia,dimensao,i,T,corpos);
			Fp(dx,r,m,distancia,i*dimensao,dimensao,T,corpos);
			ruthP(dimensao,p,dt,dx,d[j],i*dimensao,T);
		}
		if(i%intervaloGravacao==0)
		{
			escreverArquivos();
		}
	}
}

void calcularDistancia(double *v,double *distancia, int dimensao, int instante, int T,int corpos)
{
	int i,j,k;

	for(i=0;i<corpos;i++)
	{
		for(j=0;j<corpos;j++)
		{
			distancia[i*corpos+j]=0;
			for(k=0;k<dimensao;k++)
			{
				temp=(v[i*dimensao+k]-v[j*dimensao+k]);
				temp=temp*temp;
				distancia[i*corpos+j]=distancia[i*corpos+j]+temp;
			}
			distancia[i*corpos+j]=(double)sqrt((double)distancia[i*corpos+j]);
		}
	}
}

void Fr(double *dx, double *p, double *m, int inicio,int dimensao, int T)
{
	int i,j;

	// corpo
	for(i=0;i<corpos;i++)
	{
		//coordenada
		for(j=0;j<dimensao;j++)
		{
			//dx[corpo][coordenada]= dx[corpo*dimensao+coordenada]
			dx[i*dimensao+j]= p[i*dimensao+j]/m[i];
		}
	}
}

void ruthR(int dimensao, double *r,double dt,double *dx,double c,int inicio,int T)
{
	int i,j;

	// corpo
	for(i=0;i<corpos;i++)
	{
		//coordenada
		for(j=0;j<dimensao;j++)
		{
			r[i*dimensao+j]=r[i*dimensao+j]+dt*c*dx[i*dimensao+j];
		}
	}
}

void Fp(double *dx,double *r, double *m,double *distancia, int inicio, int dimensao,int T, int corpos)
{
	int i,j,k;

	// corpo
	for(i=0;i<corpos;i++)
	{
		//coordenada
		for(j=0;j<dimensao;j++)
		{
			dx[i*dimensao+j]=0;
			for(k=0;k<corpos;k++)
			{
				//se k nao for o corpo computado
				if(k!=i){
					dx[i*dimensao+j]= dx[i*dimensao+j] -(G*m[i]*m[k]*(r[i*dimensao+j]-r[k*dimensao+j]))/(double)(pow((double)distancia[i*corpos+k],3));
				}
			}
		}
	}
}

void ruthP(int dimensao, double *p,double dt,double *dx,double d, int inicio, int T)
{
	int i,j;

	// corpo
	for(i=0;i<corpos;i++)
	{
		//coordenada
		for(j=0;j<dimensao;j++)
		{
			p[i*dimensao+j]=p[i*dimensao+j]+dt*d*dx[i*dimensao+j];
		}
	}
}

void escreverArquivos(void)
{
	int j,k;
	char aux[100];

	for(k=0;k<corpos;k++)
	{
		for(j=0;j<dimensao;j++)
		{
			//arquivo txt
			//posicao
			sprintf(aux,"%0.20lf ",r[k*dimensao+j]);
			fputs(aux,f[k]);
			//momento
			sprintf(aux,"%0.20lf ",p[k*dimensao+j]);
			fputs(aux,f[k+corpos]);
		}
		fputs("\n",f[k]);
		fputs("\n",f[k+corpos]);
	}
}

void calcularTempo(void)
{
	fim=clock();
	total=fim-inicio;
	tempo=(double)total/CLOCKS_PER_SEC;
	tempo2=(int)tempo;
	horas=(tempo2/3600);
	minutos=(tempo2/60)-horas*60;
	segundos=(tempo2 % 60);
	printf("%d h %d min %d s\n",horas, minutos ,segundos);
	printf("%lf \n",tempo);
}

void fecharArquivos(void)
{
	char aux[100];

	int NFiles= corpos*2+1,i;
	sprintf(aux,"tempo de execucao= %0.3lf s",tempo);
	strcat(aux,"\n");
	fputs(aux,f[NFiles-1]);

	for(i=0;i<NFiles;i++)
	{
		fclose(f[i]);
	}
}
